Autosolitons in trapped Bose-Einstein condensates with two- and three-body inelastic 

processes 
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In this work, we consider the conditions for the existence of autosolitons, in trapped Bose-Einstein 
condensates with attractive atomic interactions. First, the variational approach is employed to es- 
timate the stationary solutions for the three-dimensional Gross-Pitaevskii equation with trap po- 
tential, linear atomic feeding from the thermal cloud and two- and three-body inelastic processes. 
Next, by using exact numerical calculations, we show that the variational approach gives reliable 
analytical results. We also discuss the possible observation of autosolitons in experiments with 7 Li. 

PACS numbers: 03.75.Fi, 32.80.Pj, 42.50.Md, 42.81.Dp 
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I. INTRODUCTION 

The existence of envelope autosolitons in one- 
dimensional (ID) case, in homogeneous nonlinear 
medium with dissipation and amplification, was revealed 
by Pereira and Stenflo JO]. They found the exact solu- 
tion for autosolitons with arbitrary growth and damp- 
ing strengths in the perturbed nonlinear Schrodinger 
equation (NLSE). Later, autosolitons were discovered 
in nonlinear fiber optics, namely in fibers with am- 
plifiers and distributed filters (the latter corresponds 
to the frequency-dependent damping in the nonlinear 
Schrodinger equation) | z ,|| and also for waves on the 
surface of deep water [ 1 1 . Correspondingly, in a two- 
dimensional (2D) homogeneous medium with amplifica- 
tion and nonlinear damping, the possibility of existence 
of a 2D analog of the Pereira-Stenflo solitons was re- 
cently shown by a variational approach Q. Autosoli- 
tons in a weakly dispersive nonlinear media, described 
by the Korteweg-deVries equation, have been studied in 
Refs. [§§. 

The autosolitons can be distinguished from ordinary 
solitons. The latter exist in conservative media and are 
originated from the balance between the nonlinear and 
dispersive effects of the wave propagation. The proper- 
ties of these generated solitons are defined by the initial 
conditions (their number, parameters like amplitudes, 
widths, etc.) ||, with the solutions characterized by 
their corresponding properties. As for to the autosoli- 
tons, they can be generated in nonconservative media 
when effects of amplification and dissipation are present. 
For the existence of autosolitons, one should add to the 
equilibrium condition between nonlinearity and disper- 
sion the requirement of a balance between amplification, 
frequency-dependent damping and nonlinear dissipation. 
In distinction from ordinary solitons, the properties of 
the autosoliton, as a rule, are fixed by the coefficients of 
the perturbed NLSE and by any initial perturbation that 
is attracted to this point (attractor in the space of coeffi- 
cients). Mathematically, the problem is described by the 
NLSE with complex parameters. If the imaginary parts 



are large, the equation is equivalent to the so-called com- 
plex Ginzburg-Landau equation. An interesting limit is 
represented by the NLSE with small complex coefficients. 

The purpose of this article is to show that the ana- 
log of autosoliton is possible in a trapped Bose-Einstein 
condensate (BEC). Recently the existence of bright and 
dark solitons in BEC has been demonstrated (see theory 
in [|]|J). Dark solitons have already been observed in 
BEC with repulsive interaction between atoms |hJ . For 
attractive interactions, bright solitons exist in ID BEC. 
It is well known that two- and three-dimensional con- 
densates with attractive interaction between atoms and 
trapping potential are unstable if the number of atoms 
exceeds the critical value (N c ). Below this value a stable 
ground state can exist, corresponding to solitary solu- 
tions Jll],[l2| . When the number of atoms exceeds N c the 
collapse occurs. At large densities the inelastic scattering 
processes involving two and three atoms come into play, 
leading to the effective nonlinear damping of the con- 
densate. The process of feeding atoms from the thermal 
cloud can be modeled as a linear amplification described 
in Ref. O] , where the relevance of the three-body inelas- 
tic processes was discussed. The statistical data obtained 
from experiments with 7 Li supports the growing and col- 
lapse picture 14 jl| . 

Numerical simulations of the 3D Gross-Pitaevskii (GP) 
equation are performed in Refs. Jl3| , ^6| . The present work 
shows that periodic oscillations occur in the condensate 
and, for particular values of atomic feeding and three- 
body dissipation parameters, stable states of the cloud 
can exist. Thus, we can expect the occurrence of analog 
autosolitons in 2D and 3D BECs. The problem is de- 
scribed by the complex Ginzburg-Landau equation with 
trapping potential in the NLSE limit with small noncon- 
servative perturbations. This equation is nonintegrable 
and, therefore, analytical solutions can only be obtained 
by considering approximate methods like the variational 
approach Jl7[ |. In the following, we first use the time- 
dependent variational approach to obtain the solutions. 
As showed previously, the time-dependent variational ap- 
proach is quite effective to study the dynamics of 3D 
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BEC in trapping potential with conservative perturba- 
tions |H . These autosoliton solutions can be considered 
as the nonlinear modes of such system like solitons for 
the integrable NLSE M|. Exact numerical simulations 
are also performed in the present work, which confirms 
that, in the present case, the variational approach is a 
convenient and reliable approximation. 



II. THE MODEL AND THE VARIATIONAL 
APPROACH 

Let us describe the dynamics of a trapped Bose- 
Einstein condensate in the framework of the Gross- 
Pitaevskii equation. The space-time dynamics of the con- 
densate wavefunction can be analyzed by the GP equa- 
tion in the mean field approximation: 

\u t + Au - (ft 2 r 2 )u + X 2 \u\ 2 u + \ 3 \u\ 4 u 

= i{ 7 - t i\u\ 2 -^\u\ 4 )u = R(u,u*), (1) 

where we use a standard simplified notation for the time 
derivative through a lower index t. ft 2 r 2 is the trap har- 
monic potential; A2 and A3 are, respectively, the two and 
three body interaction parameters, where A2(= a s ) is 
given by the s- wave atomic scattering length. 7, fi, and 
£ are positive defined coefficients related, respectively, to 
feeding, dipolar relaxation and three-body inelastic re- 
combination parameters. We drop the explicit time and 
radial dependence of the functions, unless they are nec- 
essary for clarity. 

In the present variational approach, for u = u(r, t), we 
use the Gaussian trial function [EOl 



u = A(t) exp — 



2a 2 (t) 



(2) 



where A, a, 6, and (j) are : respectively, the amplitude, 
width, chirp and linear phase. We did not include the 
center-of-mass coordinate into the ansatz, because the 
dissipative and amplifying terms have no influence on it. 

The variational approach is applied to the averaged 
Lagrangian of the conservative system 



L = J C(r,t)d 3 r, 



(3) 



where the Lagrangian density, C = C(r,t), is given by 



A 2 



A, 



C = -(u t u* - u* t u) - \Vu\ 2 + -f\u\ 4 + f\u\» - QV|h| 2 

(4) 

Substituting the trial function (||) into Eqs. (Q) and (||), 
we find the averaged Lagrangian in terms of the conden- 
sate wavefunction parameters 



L 



-A 2 a 3 



G 



3a 2 b t + 4^ + ^(1 + a 4 b 2 ) 



" 2 A 2 



4A 



3 ,4 



a/2 9V3 



A' 1 + 6ft 2 a 2 



(5) 



We formally add Cr to Eq.(Q), with the property that 
SCr/Su* = —R(u,u*), where R is the right hand side 
of Eq. (|l]) . By applying the Euler-Lagrange equations to 
£' = C + Cr, with respect to u*, we obtain 



dC 

du* 



d dC 
dt dul 



dC 

du* 



d dC 
dt duf 



R(u,u*) = 0, (6) 



that leads to Eq. ([j]) (The conjugate equation is obtained 
in a similar way). 

The corresponding variational principle is given by 

S [ L'dt = 5 [ (L+ L R )dt = 0, (7) 
Jo Jo 

where, as in Eq.(|), L R = J d 3 rC R . Taking into account 
that for a small shift Sr] of some variational parameter 77, 
we have 



df 

f(v + dv) = f(V) +St 1-q-, 



(8) 



where / = f(u,u*) = L or Lr, we obtain a system of 
equations for the variational parameters m [p0pl|] : 



8L 

dm 



d dL 
dt dr/ti 



= d 



R 



R* 



du 

dm 



(9) 



where Eq. (g) and its conjugate were used. 

The substitution of Eqs.(||) and (||) into Eq.(||), yields 
the following system of ODE's : 



d(A 2 a 3 ) 
dt 

d(A 2 a 5 ) 



2l A 2 a 3 - -^A 4 a 3 
v2 



£A 6 a 3 , 



3^3 



dt 



db 
dt 
dtp 
~dt 



AA 2 a 5 b + 2jA 2 a 5 
^- -2b 2 - 2ft 2 - 



4„5 



2^/2 

2 „,o „„t \ 2 A 2 



A^a 



6„5 



9V3 



£_A b a 



7 



■MA 2 



AX 3 A 4 
2V2a 2 9V3a 2 ' 
2 



= A,A 4 . 



8V2- ' • (10) 

TheEq. (0a) can also be obtained from the modified 
form of the conservation law for the number of atoms N, 
where N is given by 



N 



\u\ 2 d A r 



ill) 



The other equations of the system ( |l0|) can be obtained 
by using higher moments, as shown in the appendix. It 
is useful to rewrite the system using the notation x = a 2 , 
V = A 2 : 



x t = Axb 



2V2 
6yb + 2~/y 



xy 



9VS 
7/i 2 



b t = — -2b z - 2ft 2 - 



4V2 a 3x/3^ ' 
A 2 ?/ 4A 3 j/ 2 



2y/2x 9\/3a 



(12) 



This ODE system is the main result of this section. 
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III. ANALYSIS OF THE FIXED POINTS 

The autosoliton solution corresponds to the fixed 
points of the system. Recall the main properties of the 
autowave (dissipative) soliton. This solution has the form 
like for standard soliton, defined by the balance between 
amplification and nonlinear dissipative terms. As showed 
in the analysis of the one dimensional case (the Pereira- 
Stenflo soliton), the solution is fixed by the parameters 
related to the amplification and dissipation with chirped 
phase. It does not depend on the initial conditions. 

We restrict our analysis to A3 = 0. From the system 
given by Eq. (|l^), we can obtain the fixed points. In the 
next we distinguish three cases: 

1. Case fi = 0, £ ^ 0. 

In this case, 



When fi = (£ ^ 0), we recover Eqs.([L3|) and (14). Al- 
ternatively, when £ — > (fi 7^ 0) we recover Eq.(16). 

We now investigate the stability of the fixed points, in 
cases 1 and 2, by using the linear stability analysis: 

1. Case fi = 0, £ ^ 0. 

Let present the solutions of the system as x — x s \ +x\, 
y = y s i + T/i, b = b s i + b\. The linearized system for cor- 
rections is 



x u = -^-j=£,y s iXsiyi + 4x s i6i = -c 2 yi + c 3 bi, 

Vit = -872/1 - Qy s \b\ = -d 2 yi - d 3 bi, 

, 4 ( X 2 y sl 1 \ A 2 
"it = — I —ef I x l ~ ^r-^ 2/1 - 4o s i6i 



x 2 sl V 8%/2 x s J " l 2y/2x al ' 
aixi - a 2 yi - a 3 bi. 



(19) 



/3V3 7 , 1 



and the width is 



x.i = -y + Vx + fcl ' 



(13) 



(14) 



With the solutions xi, y\ and b\ having the same expo- 
nential behavior in time, given by ~ e qt , we obtain the 
characteristic equation 



q +aiq -a 2 q-a 3 =0, 



where 



(20) 



where 



Pi = . 7TT' fe l 



AV2(W + b 2 si y & + b 2 sl 



(15) 



Let us consider that p\ >> k\. Then the solution is 



x s i 



Pi A 2 j/ S i 



x 20 
a x = {d 2 + a 3 ) = y7, 

a 2 = (aic 3 + a 2 d 3 - d 2 a 3 ) 

32 2 ( \ 2 y s ix sl _ _\ 3\ 2 y sl 

~ 3 7 x 2 sl \ 8V2 J V2x sl 

a 3 = ai(d 3 c 2 + d 2 c 3 ) = I 1 ] . ( 21 ) 



2. Case £ = 0, fi ^ 0. 

The fixed points are 



The roots with Re(q) > correspond to the unstable 
solutions. 



2%/2 7 7 

Vs2 = , b s2 = --. 

/i 4 



(16) 



and the width is defined by the same expression as before 
Eq.Qwith y s2 ,b s2 . 

3. Case fi ^ 0, £ ^ 0. 

The fixed points are: 



Vs3 



3V3fi / 27/i 2 3-y/3 7 
V 32^ + ^' 



, M 7 



where 



A2y. 



«3 



1 



4V2(fi 2 +62 3 )' 



ft 2 + fe 2 



(17) 



(18) 
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FIG. 1. Gaussian variational analysis of stability of the 
fixed points for the Gross-Pitaevskii equation including feed- 
ing (7) and three-body losses (£). 
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Without loss of generality, in our dimensionless NLSE 
we can scale the parameters as X 2 = 1 and fi = 1 [fl3|| . 
The diagram of stability, according to the solutions of 
Eq. (|2(i|), is presented in Fig. 1. The diagram clearly 
shows that, when 7 > 1.84£, the system is unstable. If 
7 >> £, the system enters the collapsing process that has 
been shown in Ref. Jl6| to be chaotic. If 7 is decreased 
(or £ increased) the system will eventually achieve a sta- 
ble region where the formation of autosoliton is possible. 



With the same scaling used in case 1 (A2 = 1 and 
ft = 1), and with the above solutions of Eq.(^3|), we 
obtain the diagram of stability, shown in Fig. 2. The 
diagram clearly shows that, when 7 > 0.53/i, the sys- 
tem is unstable. In analogy with the previous case, if 7 
is decreased (or fi increased), the system will eventually 
achieve a stable region where the formation of autosoliton 
is possible. 



2. Case £ = 0,/x ^ 0. 

Analogously, as in case 1, we present the solutions of 
the system as x — x s2 +x 2 , y = y S 2 + y2, and b = b s2 + b 2 . 
The linearized system for corrections is given by 



X2t 



V2t 



t-LX s2 



2\/2 



2/2 + 4:X s2 b2, 



:1V2 



12V27, 



ht = -5- ( -P-Xs2 - 1 ) X 2 



S2 



4// 



2\/2: 



X S 2 



-2/2-46 s2 6 2 . (22) 



The system has the same form as in Eq. (|19|). Corre- 
spondingly, the characteristic equation is given by 



q 6 + fcq 2 - (3 2 q ~ fa = 0, 



(23) 



where 



01 
02 

03 



2 " 

7 2 
-7 

2 ' 



16 / A 2 7x s 



's2 



4// 



- 1 



6A27 

X S 2H ' 



327 ( \2rXs 
17 



(24) 
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IV. NUMERICAL SIMULATIONS AND 
DISCUSSIONS 

We did a series of time-dependent simulations of the 
system within the Gaussian variational approach, using 
Eq. (|l2"|) , and also by performing exact numerical calcu- 
lations with Eq. (jl|) . In our numerical calculations, we 
have used the finite-difference Crank- Nicolson algorithm. 
The exact initial wave functions were used following the 
prescription given in |2^j . In the next, we present sim- 
ulations in a range of parameters that lead to long-time 
stable autosolitons. We obtain results with autosolitonic 
solutions for a class of parameters that are near the re- 
alistic ones, as indicated by 7 Li experiments in BEC. 
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FIG. 2. Gaussian variational analysis of stability of the 
fixed points for the Gross-Pitaevskii equation including feed- 
ing (7) and two-body losses(/i). 
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FIG. 3. Evolution of the number of atoms N in the 
Gross-Pitaevskii equation with feeding parameter 7 = 1CP 3 
and three-body dissipation parameter £ = (/x = 0). The 

results are represented in dashed lines for the variational ap- 
proach, and in solid lines for the exact numerical calculations. 
Cases I, II and III correspond, respectively, to the initial con- 
ditions N(t = 0)/iV c =0.1, 0.2 and 0.75, where iV c is the max- 
imum critical number for stability, t is given in units of the 
inverse of the trap frequency Q. 

In Fig. 3, for 7 = 1CP 3 and £ = 10~ 3 , we show the 
time evolution of the number of atoms, in terms of the 
maximum critical number for stability, N c . The forma- 
tion of the autosoliton is demonstrated either by Gaus- 
sian variational approach or by exact numerical calcu- 
lations. There is a remarkable agreement between both 
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approaches. Note that the number of atoms does not 
depend on the initial conditions, but is related to the 
equilibration of the feeding and dissipation. The vari- 
ational approach results give a little higher number of 
atoms than the exact calculations. 
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FIG. 4. Evolution of the mean square radius, Eq. (|25|), 
in the Gross Pitaevskii equation. The same parameters and 
conventions given in Fig. 3 were used. 

Corresponding to the results in Fig. 3, we show in Fig. 
4 the results for the time evolution of the mean square 
radius, where 



(r 2 ) 



r l \v,\ A (Tr. 



(25) 



In this case, the variational approach results are a bit 
lower than the ones obtained by exact calculations. 

In analogy with the case that fi = 0, represented in 
Figs. 3 and 4, we also present results for the case that the 
three-body dissipation parameter (£) is zero. The results 
obtained for the time evolution of the number of atoms 
and the mean-square-radius are, respectively, shown in 
Figs. 5 and 6, for 7 = 5 x 1CT 5 and p = 1CT 4 . The 
formation of the autosoliton is also demonstrated either 
by Gaussian variational approach or by exact numerical 
calculations. The remarkable agreement between both 
approaches, already observed in case that /1 = (Figs. 
3 and 4), also apply to this case that £ = 0, with the 
number of atoms not depending on the initial conditions. 

In 7 Li experiment the feeding parameter can be in- 
directly inferred from measurements done by the Rice 
Group jl4| and will correspond to an average rate of 
about 600 atoms/s p3|l . The two- and three-body losses 
were also measured ]24j] and estimated pq| , giving atom 
loss rates of about 2 x 10 _14 cm 3 /s and ~ 10~ 28 cm 6 /s. 
These rates were measured for non-condensed atoms. For 
condensed atoms they must be divided by factors of 2! 
and 3!, respectively With a scaled equation, such 

that A2 = 1 and Q = 1, we have for condensed atoms 
7 - 10" 3 , n ~ 10~ 4 and £ - 10~ 6 . Considering this 



magnitudes, the autosoliton can possibly be observed ex- 
perimentally in 7 Li, either by results with decreasing 7 
due to losses, or by increasing the dissipation rate due to 
other mechanism as the Feshbach resonances P7| . In case 
of diminishing 7, the autosoliton formation is more likely 
determined by dipolar relaxation rather than three-body 
recombination losses. 
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FIG. 5. Evolution of the number of atoms in the 
Gross-Pitaevskii equation with feeding 7 = 5 x 1CP J and 
two-body dissipation parameter /i — 1CP 4 (£ = 0). The initial 
conditions and conventions are the same as in Fig. 3. 
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FIG. 6. Evolution of the mean square radius in the Gross 
Pitaevskii equation. The parameters and conventions are the 
same as in Fig. 5. 



V. CONCLUSION 

In this work we study the possibility of existence 
of autosolitons in trapped 3D BEC in the presence of 
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two- and three-body inelastic processes: dipolar relax- 
ation and three-body recombination. Using the time- 
dependent variational approach for nonconservative 3D 
Gross-Pitaevskii equation, we derive expressions for the 
parameters of autosoliton and check their stability. The 
results obtained by using the present time-dependent 
Gaussian variational approach, in the NLSE with atomic 
feeding and nonlinear dissipative terms, show a remark- 
able agreement with exact numerical calculations, when 
the parameters were such that stabilization was achieved. 
We should note that, stabilization can also be numeri- 
cally observed in the model given in Ref . Jl3| ] , where 7 is 
decreased in time, but their parameters are far from the 
realistic ones, and dipolar relaxation was neglected. Such 
stabilization verified in Ref. jl3) has not been recognized 
as a manifestation of autosoliton, a characterization that 
we are pointing out in the present work. 

We have shown that the transition from unstable (col- 
lapsing) to stable point (autosoliton) solely depends on 
the magnitude of the parameters. Also the present work 
include non negligible two-body dissipative effects that 
model the dipolar relaxation losses, and that can be asso- 
ciated with values measured in ultracold 7 Li . In case 
of decreasing 7 due to collapsing and losses, the autosoli- 
ton is more likely to be formed first due to dipolar relax- 
ation rather than by three-body recombination processes. 
These results can be relevant in current experiments with 
attractive scattering length and possibly display new phe- 
nomenon of Pereira-Stenflo type autosoliton formation in 
Bose condensates. We believe that such phenomenon is 
already occurring in the long time behavior in the ac- 
tual experiments with 7 Li jl4j (fi ~ 2nx 140Hz), since 
for longer times (~60s) the maximum number of atoms 
(N c ~ 1300 atoms) is considerably reduced, as expected 
in our simulations. We hope that experiments with di- 
rect observation of the evolution of the condensate can 
clarify this matter. 



APPENDIX 

The first equation of the system ( ]To| ) can be obtained 
by calculating the rate of change of the number of atoms 
as 



dN _ d 
~dt ~ dt 



= Jl / rf 3 , ,2 



{u* t u + u*ut)d 3 r. 



(26) 



Taking u t from Eq. (|l|) and substituting in ( |26| ) we obtain 
the modified form of the conservation law for the number 



of atoms 1 29 



^ = 2 7 iV - 2^i J \u\ 4 d 3 r -2$J \u\ 6 d 3 r. (27) 

Substituting into this equation the Gaussian trial func- 
tion (|^) we obtain the first equation of the system given 
in Eq. ©. 



We have derived the Eqs. (p. Op) and (|10|c) by the mo- 
ments method. Eq. (|Io|b) can be derived by calculating 



d(r 2 ) 



dt 



1 * 2 * 2 \ 

(u t r u-\-u r u t )a r. 



(28) 



Substituting the Ut from expression (|l]) and applying the 
commutation rules, we get |2q] 



d(r 2 ) r r 

' = 4Im / u*f- (Vu)tfV + 2-y \u\ 2 r 2 d 3 r 

^/ N W r -2f/ HW , (29, 



dt 



The substitution of the Gaussian ansatz (Eq. 2j) in both 
sides of this expression results directly in Eq.(lGb). 

Equation (10:) can be derived analogously by showing 
that 

= -4f2 2 Im J u*r- (X7u)d 3 r 

+2Im J (Vu*) • V(V lU )d 3 r 

+2Re J (Vtt*) • V(V 2 u)d 3 r (30) 

where V\ — — A 2 |u| 2 — A 3 |-u| 4 and V 2 = 7 — ^M 2 — £|w| 4 - 
Substituting the Gaussian ansatz (Eq. ||) in both sides 
of Eq. ( [BO] ) and using results (|i"o|a) and (|l0|b) there is 
overall cancellation of the feeding and dissipative terms 
and we finally get Eq. (nT 
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